The R code from this presentation is here: https://raw.githubusercontent.com/IvanVoronin/TwinAnalysis/master/slides/intro/intro.R
# Making univariate toy data
library(dplyr,
quietly = TRUE)
mz_data <- MASS::mvrnorm(
n = 100,
mu = c(0, 0),
Sigma = matrix(c(1.0, 0.8,
0.8, 1.0),
nrow = 2)) %>%
as.data.frame() %>%
setNames(c('X1', 'X2'))
dz_data <- MASS::mvrnorm(
n = 100,
mu = c(0, 0),
Sigma = matrix(c(1.0, 0.5,
0.5, 1.0),
nrow = 2)) %>%
as.data.frame() %>%
setNames(c('X1', 'X2'))
data <- rbind(data.frame(zyg = 'MZ', mz_data),
data.frame(zyg = 'DZ', dz_data)) %>%
mutate(
zyg = factor(zyg, levels = c('MZ', 'DZ'))
)
head(data)
## zyg X1 X2
## 1 MZ 1.6139277 1.7912405
## 2 MZ -0.7454432 -0.4165514
## 3 MZ -1.2403718 -0.6782260
## 4 MZ -0.9808771 -1.0968942
## 5 MZ -0.3770839 -0.2187076
## 6 MZ -0.6517744 -0.4693972# Installing the packages
# install.packages('devtools')
# library(devtools)
# install_github('ivanvoronin/mlth.data.frame')
# install_github('ivanvoronin/TwinAnalysis')
# Loading the packages
library(OpenMx)
library(TwinAnalysis)
# Descriptive statistics
get_univ_descriptives(data, zyg = 'zyg', vars = 'X')
## $X
## Twin1-------------------- Twin2-------------------- Cor--------------------------
## n M SD n M SD r lowerCI upperCI
## MZ 100 -0.06778361 0.9888608 100 -0.075597225 1.022780 0.7881512 0.7001161 0.8525833
## DZ 100 0.06555748 1.0586309 100 -0.008591523 1.234908 0.5030154 0.3402106 0.6365423# Univariate ACE model
model <- univariate_ace(data, zyg = 'zyg', ph = 'X')
model <- mxRun(model)
summary(model)
## Summary of ACE
##
## free parameters:
## name matrix row col Estimate Std.Error A
## 1 ACE.mean[1,1] mean 1 1 -0.01625766 0.06929624
## 2 ACE.a[1,1] a X X 0.91070049 0.09748841
## 3 ACE.c[1,1] c X X 0.37057316 0.22770863
## 4 ACE.e[1,1] e X X 0.46506707 0.03335832
##
## Model Statistics:
## | Parameters | Degrees of Freedom
## Model: 4 396
## Saturated: NA NA
## Independence: NA NA
## | Fit (-2lnL units)
## Model: 1067.62
## Saturated: NA
## Independence: NA
## Number of observations/statistics: 200/400
##
## Information Criteria:
## | df Penalty | Parameters Penalty
## AIC: 275.620 1075.620
## BIC: -1030.514 1088.813
## | Sample-Size Adjusted
## AIC: 1075.825
## BIC: 1076.141
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 04:22:46
## Wall clock time: 0.07753515 secs
## optimizer: SLSQP
## OpenMx version number: 2.18.1
## Need help? See help(mxSummary)# Constrained models
# This is a list of MxModels
nested_models <- def_nested_models(
model,
AE = mxConstraint(c[1, 1] == 0),
CE = mxConstraint(a[1, 1] == 0),
E = list(mxConstraint(c[1, 1] == 0),
mxConstraint(a[1, 1] == 0)),
run = TRUE
)
fit_stats(model, nested_models = nested_models)
## Warning in computeFitStatistics(likelihood, DoF, chi,
## chiDoF, retval[["numObs"]], : Your model may be mis-
## specified (and fit worse than an independence model),
## or you may be using the wrong independence model, see ?
## mxRefModels
## model base ep minus2LL df AIC BIC
## 1 ACE Saturated 4 1067.620 396 275.6200 -1030.5137
## 2 AE ACE 4 1068.249 397 274.2487 -1035.1833
## 3 CE ACE 4 1095.945 397 301.9449 -1007.4871
## 4 E ACE 4 1194.151 398 398.1515 -914.5788
## CFI TLI RMSEA diffLL diffdf
## 1 0.992978093 0.9976594 0.02696322 6.8724180 6
## 2 0.995966781 0.9988477 0.01891893 0.6286784 1
## 3 0.773045631 0.9351559 0.14191883 28.3249192 1
## 4 -0.009349393 0.7476627 0.27995972 126.5314918 2
## p
## 1 3.328073e-01
## 2 4.278405e-01
## 3 1.025671e-07
## 4 3.342225e-28# Making bivariate toy data
mz_data2 <- MASS::mvrnorm(
n = 100,
mu = c(0, 0, 0, 0),
Sigma = matrix(
c(1.0, 0.7, 0.6, 0.5,
0.7, 1.0, 0.5, 0.8,
0.6, 0.5, 1.0, 0.7,
0.5, 0.8, 0.7, 1.0),
nrow = 4)) %>%
as.data.frame() %>%
setNames(c('X1', 'Y1', 'X2', 'Y2'))
dz_data2 <- MASS::mvrnorm(
n = 100,
mu = c(0, 0, 0, 0),
Sigma = matrix(
c(1.0, 0.7, 0.4, 0.2,
0.7, 1.0, 0.2, 0.4,
0.4, 0.2, 1.0, 0.7,
0.2, 0.4, 0.7, 1.0),
nrow = 4)) %>%
as.data.frame() %>%
setNames(c('X1', 'Y1', 'X2', 'Y2'))
data2 <- rbind(data.frame(zyg = 'MZ', mz_data2),
data.frame(zyg = 'DZ', dz_data2)) %>%
mutate(
zyg = factor(zyg, levels = c('MZ', 'DZ'))
)
head(data2)
## zyg X1 Y1 X2 Y2
## 1 MZ -1.0800889 -1.0275203 -0.1378594 -0.82322342
## 2 MZ 0.8290774 0.9620859 0.9405737 0.88065324
## 3 MZ 0.5680403 0.1000585 1.2833039 0.41037703
## 4 MZ 0.2603521 1.2183755 -0.4669039 0.90408497
## 5 MZ -0.1713043 -0.4440812 0.7967614 0.04784216
## 6 MZ -0.2450042 1.0047931 -0.7089274 0.50509900# Descriptive statistics
get_univ_descriptives(data2,
zyg = 'zyg',
vars = c('X', 'Y'))
## $X
## Twin1------------------ Twin2-------------------- Cor--------------------------
## n M SD n M SD r lowerCI upperCI
## MZ 100 -0.1359043 1.003993 100 -0.06795439 0.9190689 0.5930486 0.4489236 0.7070999
## DZ 100 -0.1636363 1.017908 100 -0.13654309 0.9801215 0.5075549 0.3455896 0.6401541
##
## $Y
## Twin1-------------------- Twin2------------------- Cor--------------------------
## n M SD n M SD r lowerCI upperCI
## MZ 100 -0.02955259 0.9583681 100 -0.1114836 0.9391590 0.7579110 0.6597008 0.8306695
## DZ 100 -0.21072036 1.1037023 100 -0.1693055 0.9790516 0.5489379 0.3951225 0.6728125# Cross-twin cross-trait correlations
get_cross_trait_cors(data2,
zyg = 'zyg',
vars = c('X', 'Y'))
## zyg: MZ
## X1 Y1 X2 Y2
## X1 1.0000000 0.7121176 0.5930486 0.4323049
## Y1 0.7121176 1.0000000 0.5325547 0.7579110
## X2 0.5930486 0.5325547 1.0000000 0.7156957
## Y2 0.4323049 0.7579110 0.7156957 1.0000000
## --------------------------------------------
## zyg: DZ
## X1 Y1 X2 Y2
## X1 1.0000000 0.7162098 0.5075549 0.3259307
## Y1 0.7162098 1.0000000 0.3067146 0.5489379
## X2 0.5075549 0.3067146 1.0000000 0.6645038
## Y2 0.3259307 0.5489379 0.6645038 1.0000000# Bivariate ACE model
model2 <- multivariate_ace(
data = data2,
zyg = 'zyg',
ph = c('X', 'Y')
)
model2 <- mxRun(model2)
summary(model2)
## Summary of ACE
##
## free parameters:
## name matrix row col Estimate Std.Error A
## 1 ACE.mean[1,1] mean 1 1 -1.291566e-01 0.06079859
## 2 ACE.mean[1,2] mean 1 2 -1.348547e-01 0.06403822
## 3 ACE.a[1,1] a X X 1.897923e-06 0.28839284
## 4 ACE.a[2,1] a Y X 4.925751e-01 0.12721120
## 5 ACE.a[2,2] a Y Y 7.182836e-01 0.09900747
## 6 ACE.c[1,1] c X X 5.335468e-01 0.04600191
## 7 ACE.c[2,1] c Y X 2.438873e-01 0.18791282
## 8 ACE.c[2,2] c Y Y 5.127795e-01 0.13539231
## 9 ACE.e[1,1] e X X 4.191573e-01 0.02227449
## 10 ACE.e[2,1] e Y X 4.407855e-01 0.04928814
## 11 ACE.e[2,2] e Y Y 4.701544e-01 0.03325551
##
## Model Statistics:
## | Parameters | Degrees of Freedom
## Model: 11 789
## Saturated: NA NA
## Independence: NA NA
## | Fit (-2lnL units)
## Model: 1752.1
## Saturated: NA
## Independence: NA
## Number of observations/statistics: 200/800
##
## Information Criteria:
## | df Penalty | Parameters Penalty
## AIC: 174.0997 1774.100
## BIC: -2428.2727 1810.381
## | Sample-Size Adjusted
## AIC: 1775.504
## BIC: 1775.532
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 04:22:46
## Wall clock time: 0.08861446 secs
## optimizer: SLSQP
## OpenMx version number: 2.18.1
## Need help? See help(mxSummary)# Output tables from the bivariate ACE
get_output_tables(model2)
## $Variance_components
## A(%) C(%) E(%)
## X 0.2535936 0.3597038 0.3867026
## Y 0.5159729 0.2629640 0.2210630
##
## $Covariation_components
## Total A C E A(%)
## X <-> X 0.9567681 0.2426303 0.3441531 0.3699847 0.2535936
## X <-> Y 0.6861063 0.3538086 0.1250604 0.2072373 0.5156761
## Y <-> Y 0.9999194 0.5159313 0.2629428 0.2210452 0.5159729
## C(%) E(%)
## X <-> X 0.3597038 0.3867026
## X <-> Y 0.1822756 0.3020483
## Y <-> Y 0.2629640 0.2210630
##
## $All_correlations
## r_A r_C r_E Total
## X <-> X 1 1.0000000 1.0000000 1.0000000
## X <-> Y 1 0.4157318 0.7246619 0.7014643
## Y <-> Y 1 1.0000000 1.0000000 1.0000000A model (mxModel):
# Twin models in OpenMx
library(OpenMx)
# Univariate ACE model
model <- mxModel(
name = 'ACE',
# Means
mxMatrix(type = 'Full',
nrow = 1, ncol = 1,
free = TRUE,
name = 'mean'),
mxAlgebra(cbind(mean, mean),
name = 'expMeans'),
# Variance components
mxMatrix(type = 'Lower',
nrow = 1, ncol = 1,
free = TRUE,
name = 'a'),
mxMatrix(type = 'Lower',
nrow = 1, ncol = 1,
free = TRUE,
name = 'c'),
mxMatrix(type = 'Lower',
nrow = 1, ncol = 1,
free = TRUE,
name = 'e'),
mxAlgebra(t(a) %*% a, name = 'A'),
mxAlgebra(t(c) %*% c, name = 'C'),
mxAlgebra(t(e) %*% e, name = 'E'),
# Expected covariance
mxAlgebra(rbind(cbind(A + C + E, A + C ),
cbind(A + C, A + C + E)),
name = 'expCovMZ'),
mxAlgebra(rbind(cbind(A + C + E, .5%x%A + C),
cbind(.5%x%A + C, A + C + E)),
name = 'expCovDZ'),
# MZ submodel
mxModel(name = 'MZ',
mxData(observed = mz_data,
type = 'raw'),
mxExpectationNormal(covariance = 'ACE.expCovMZ',
means = 'ACE.expMeans',
dimnames = c('X1', 'X2')),
mxFitFunctionML()),
# DZ submodel
mxModel(name = 'DZ',
mxData(observed = dz_data,
type = 'raw'),
mxExpectationNormal(covariance = 'ACE.expCovDZ',
means = 'ACE.expMeans',
dimnames = c('X1', 'X2')),
mxFitFunctionML()),
mxFitFunctionMultigroup(c('MZ','DZ')),
# Output tables
mxAlgebra(A + C + E, name = 'V'),
mxAlgebra(cbind(V, A, C, E),
dimnames = list('X', c('Total', 'A', 'C', 'E')),
name = 'Raw_variance'),
mxAlgebra(cbind(A / V, C / V, E / V),
dimnames = list('X', c('A(%)', 'C(%)', 'E(%)')),
name = 'Variance_components')
)
model <- mxRun(model)
summary(model)
## Summary of ACE
##
## free parameters:
## name matrix row col Estimate Std.Error A
## 1 ACE.mean[1,1] mean 1 1 -0.01625702 0.06929625
## 2 ACE.a[1,1] a 1 1 0.91069985 0.09748699
## 3 ACE.c[1,1] c 1 1 0.37057422 0.22770404
## 4 ACE.e[1,1] e 1 1 0.46506705 0.03335826
##
## Model Statistics:
## | Parameters | Degrees of Freedom
## Model: 4 396
## Saturated: NA NA
## Independence: NA NA
## | Fit (-2lnL units)
## Model: 1067.62
## Saturated: NA
## Independence: NA
## Number of observations/statistics: 200/400
##
## Information Criteria:
## | df Penalty | Parameters Penalty
## AIC: 275.620 1075.620
## BIC: -1030.514 1088.813
## | Sample-Size Adjusted
## AIC: 1075.825
## BIC: 1076.141
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 04:22:47
## Wall clock time: 0.04803014 secs
## optimizer: SLSQP
## OpenMx version number: 2.18.1
## Need help? See help(mxSummary)# Bivariate ACE model
model2 <- mxModel(
name = 'ACE',
# Means
mxMatrix(type = 'Full',
nrow = 1, ncol = 2,
free = TRUE,
name = 'mean'),
mxAlgebra(cbind(mean, mean),
name = 'expMeans'),
# Variance components
mxMatrix(type = 'Lower',
nrow = 2, ncol = 2,
free = TRUE,
name = 'a'),
mxMatrix(type = 'Lower',
nrow = 2, ncol = 2,
free = TRUE,
name = 'c'),
mxMatrix(type = 'Lower',
nrow = 2, ncol = 2,
free = TRUE,
name = 'e'),
mxAlgebra(t(a) %*% a, name = 'A'),
mxAlgebra(t(c) %*% c, name = 'C'),
mxAlgebra(t(e) %*% e, name = 'E'),
# Covariance
mxAlgebra(rbind(cbind(A + C + E, A + C ),
cbind(A + C, A + C + E)),
name = 'expCovMZ'),
mxAlgebra(rbind(cbind(A + C + E, .5%x%A + C),
cbind(.5%x%A + C, A + C + E)),
name = 'expCovDZ'),
# MZ submodel
mxModel(name = 'MZ',
mxData(observed = mz_data2,
type = 'raw'),
mxExpectationNormal(covariance = 'ACE.expCovMZ',
means = 'ACE.expMeans',
dimnames = c('X1', 'Y1', 'X2', 'Y2')),
mxFitFunctionML()),
# DZ submodel
mxModel(name = 'DZ',
mxData(observed = dz_data2,
type = 'raw'),
mxExpectationNormal(covariance = 'ACE.expCovDZ',
means = 'ACE.expMeans',
dimnames = c('X1', 'Y1', 'X2', 'Y2')),
mxFitFunctionML()),
mxFitFunctionMultigroup(c('MZ','DZ')),
# Output tables
mxAlgebra(A + C + E, name = 'V'),
mxAlgebra(abs(A) + abs(C) + abs(E), name = 'Vabs'),
mxAlgebra(abs(A) / Vabs, name = 'Aperc'),
mxAlgebra(abs(C) / Vabs, name = 'Cperc'),
mxAlgebra(abs(E) / Vabs, name = 'Eperc'),
mxAlgebra(cbind(diag2vec(Aperc),
diag2vec(Cperc),
diag2vec(Eperc)),
dimnames = list(c('X', 'Y'), c('A(%)', 'C(%)', 'E(%)')),
name = 'Variance_components'),
mxAlgebra(cbind(vech(V),
vech(A), vech(C), vech(E),
vech(Aperc), vech(Cperc), vech(Eperc)),
dimnames = list(c('X <-> X', 'X <-> Y', 'Y <-> Y'),
c('Total', 'A', 'C', 'E', 'A(%)', 'C(%)', 'E(%)')),
name = 'Covariation_components'),
mxMatrix(type = "Iden",
nrow = 2, ncol = 2,
name = "I"),
mxAlgebra(sqrt(I * V), name = "SD_total"),
mxAlgebra(sqrt(I * A), name = 'SD_A'),
mxAlgebra(sqrt(I * C), name = 'SD_C'),
mxAlgebra(sqrt(I * E), name = 'SD_E'),
mxAlgebra(solve(SD_total) %&% V, name='r_total'),
mxAlgebra(solve(SD_A) %&% A, name = 'r_A'),
mxAlgebra(solve(SD_C) %&% C, name = 'r_C'),
mxAlgebra(solve(SD_E) %&% E, name = 'r_E'),
mxAlgebra(cbind(vech(r_A), vech(r_C), vech(r_E),
vech(r_total)),
dimnames = list(c('X <-> X', 'X <-> Y', 'Y <-> Y'),
c('r_A', 'r_C', 'r_E', 'Total')),
name='All_correlations')
)
model2 <- mxRun(model2)
summary(model2)
## Summary of ACE
##
## free parameters:
## name matrix row col Estimate Std.Error A
## 1 ACE.mean[1,1] mean 1 1 -1.291565e-01 0.06079855
## 2 ACE.mean[1,2] mean 1 2 -1.348544e-01 0.06403800
## 3 ACE.a[1,1] a 1 1 2.990996e-07 0.28843910
## 4 ACE.a[2,1] a 2 1 4.925745e-01 0.12723097
## 5 ACE.a[2,2] a 2 2 7.182822e-01 0.09902277
## 6 ACE.c[1,1] c 1 1 5.335470e-01 0.04600363
## 7 ACE.c[2,1] c 2 1 2.438879e-01 0.18794570
## 8 ACE.c[2,2] c 2 2 5.127802e-01 0.13541217
## 9 ACE.e[1,1] e 1 1 4.191574e-01 0.02227458
## 10 ACE.e[2,1] e 2 1 4.407853e-01 0.04928897
## 11 ACE.e[2,2] e 2 2 4.701545e-01 0.03325635
##
## Model Statistics:
## | Parameters | Degrees of Freedom
## Model: 11 789
## Saturated: NA NA
## Independence: NA NA
## | Fit (-2lnL units)
## Model: 1752.1
## Saturated: NA
## Independence: NA
## Number of observations/statistics: 200/800
##
## Information Criteria:
## | df Penalty | Parameters Penalty
## AIC: 174.0997 1774.100
## BIC: -2428.2727 1810.381
## | Sample-Size Adjusted
## AIC: 1775.504
## BIC: 1775.532
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 04:22:47
## Wall clock time: 0.09123135 secs
## optimizer: SLSQP
## OpenMx version number: 2.18.1
## Need help? See help(mxSummary)mxEval(Covariation_components, model2)
## Total A C E A(%)
## X <-> X 0.9567679 0.2426296 0.3441537 0.3699846 0.2535930
## X <-> Y 0.6861056 0.3538075 0.1250609 0.2072372 0.5156750
## Y <-> Y 0.9999182 0.5159293 0.2629435 0.2210453 0.5159716
## C(%) E(%)
## X <-> X 0.3597045 0.3867025
## X <-> Y 0.1822765 0.3020486
## Y <-> Y 0.2629651 0.2210634library(mlth.data.frame)
# mlth.data.frame - multiheaded data.frame
model2 <- multivariate_ace(
data = data2,
zyg = 'zyg',
ph = c('X', 'Y')
)
model2 <- def_ci(model2, ci = 'Variance_components')
model2 <- mxRun(model2, intervals = TRUE)
get_output_tables(model2)
## $Variance_components
## A(%)-------------------------- C(%)-------------------------- E(%)-------------------------
## est lCI uCI est lCI uCI est lCI uCI
## X 0.2535936 0.05284689 0.5694850 0.3597038 0.070988727 0.552154 0.3867026 0.2972783 0.4942010
## Y 0.5159729 0.26094191 0.7996993 0.2629640 0.001319719 0.492202 0.2210630 0.1623230 0.3034753
##
## $Covariation_components
## Total A C E A(%)
## X <-> X 0.9567681 0.2426303 0.3441531 0.3699847 0.2535936
## X <-> Y 0.6861063 0.3538086 0.1250604 0.2072373 0.5156761
## Y <-> Y 0.9999194 0.5159313 0.2629428 0.2210452 0.5159729
## C(%) E(%)
## X <-> X 0.3597038 0.3867026
## X <-> Y 0.1822756 0.3020483
## Y <-> Y 0.2629640 0.2210630
##
## $All_correlations
## r_A r_C r_E Total
## X <-> X 1 1.0000000 1.0000000 1.0000000
## X <-> Y 1 0.4157318 0.7246619 0.7014643
## Y <-> Y 1 1.0000000 1.0000000 1.0000000# Formatting the table for html or latex report
# Powered by knitr and kableExtra
library(kableExtra)
M %>%
behead() %>%
kable(digits = 3,
align = 'c') %>%
add_complex_header_above(M)|
A(%)
|
C(%)
|
E(%)
|
|||||||
|---|---|---|---|---|---|---|---|---|---|
| est | lCI | uCI | est | lCI | uCI | est | lCI | uCI | |
| X | 0.254 | 0.053 | 0.569 | 0.360 | 0.071 | 0.552 | 0.387 | 0.297 | 0.494 |
| Y | 0.516 | 0.261 | 0.800 | 0.263 | 0.001 | 0.492 | 0.221 | 0.162 | 0.303 |
# Saving the table for the output
M %>%
register_output(
name = 'Table 1',
caption = 'Table 1: Variance components')
## A(%)-------------------------- C(%)-------------------------- E(%)-------------------------
## est lCI uCI est lCI uCI est lCI uCI
## X 0.2535936 0.05284689 0.5694850 0.3597038 0.070988727 0.552154 0.3867026 0.2972783 0.4942010
## Y 0.5159729 0.26094191 0.7996993 0.2629640 0.001319719 0.492202 0.2210630 0.1623230 0.3034753